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ABSTRACT. Bering Glacier, Alaska, USA, has a ~20 year surge cycle, with its most recent surge reaching 
the terminus in 201 1 . To study this most recent activity a time series of ice velocity maps was produced by 
applying optical feature-tracking methods to Landsat-7 ETM+ imagery spanning 2001-11. The velocity 
maps show a yearly increase in ice surface velocity associated with the down-glacier movement of a surge 
front. In 2008/09 the maximum ice surface velocity was 1.5 ±0.01 7 km a -1 in the mid-ablation zone, 
which decreased to 1.2 ±0.01 5 km a -1 in 2009/10 in the lower ablation zone, and then increased to 
nearly 4.4 ±0.03 km a 1 in summer 2011 when the surge front reached the glacier terminus. The surge 
front propagated down-glacier as a kinematic wave at an average rate of 4.4 ±2.0 km a -1 between 
September 2002 and April 2009, then accelerated to 13.9 ±2.0 km a' 1 as it entered the piedmont lobe 
between April 2009 and September 201 0. The wave seems to have initiated near the confluence of Bering 
Glacier and Bagley Ice Valley as early as 2001, and the surge was triggered in 2008 further down-glacier 
in the mid-ablation zone after the wave passed an ice reservoir area. 


INTRODUCTION 
Bering Glacier surge history 

During the 20th century Bering Glacier, Alaska, USA 
(Fig. 1), surged five times, in ~1900, ~1920, ~1 938-40, 
1957-67 and 1993-95: approximately every 20 years 
(Molnia and Post, 2010). The 1993-95 surge was keenly 
studied using aerial photography (Lingle and others, 1993; 
Herzfeld and Mayer, 1997) and synthetic aperture radar 
(SAR) data (Fatland and Lingle, 1998; Lingle and Fatland, 
2003; Roush and others, 2003). The surge started well 
below the equilibrium-line altitude (ELA) during winter, and 
was associated with high englacial water pressures, as 
evidenced by artesian, clear, blue water in fresh crevasses 
(Lingle and others, 1993), and by interferometric detection 
of vertical displacements in Bagley Ice Valley indicating 
localized subglacial hydraulic jacking (Fatland and Lingle, 
1998; Lingle and Fatland, 2003). Analysis of SAR images 
showed the surge front propagated down-glacier at a rate of 
75-1 12 md -1 (27-40 km a -1 ) in a wave-like fashion, while 
the ice velocity was considerably less, 10-20md~' (3.7— 
7.3 km a -1 ) (Roush and others, 2003). The surge also 
propagated up-glacier at 200-500 md -1 (73-1 82 km a' 1 ) 
into Bagley Ice Valley, significantly faster than the rate at 
which it moved down-glacier. Each stage of the surge 
ended with an outburst of turbid water into proglacial Vitus 
Lake (August 1994 and September 1995) (Fatland and 
Lingle, 1998). 

Optical satellite imagery of the 1993-95 surge is not 
available due to lack of acquisitions and cloud cover, and 
SAR interferometry is only capable of producing glacier 
velocity fields over short time-spans (days to weeks), not 
over an entire year, so the pre-surge build-up and climax of 
the surge has never been studied on an annual basis. 
During the 2000s, though, Bering Glacier began a new 
surge phase and there is sufficient optical satellite imagery 
to study the surge through the entire process in the 
ablation zone. 


Previous observations of kinematic waves on glaciers 

Traveling waves on glaciers have been reported since the 
1 890s (Sharp, 1 954), with observations on the Mer de Glace, 
France, in 1891-95 (Val lot, 1900), on glaciers in Yakutat 
Bay, Alaska, in the early 1900s (Tarr and Martin, 1914), on 
Hintereisferner, Austria, in 1905 (Blilmcke and Finster- 
walder, 1905), on Black Rapids Glacier, Alaska, in 1 936— 
37 (Hance, 1937), and on Nisqually Glacier, Washington, 
USA, in the 1950s and 1960s Oohnson, 1953; Harrison, 
1956; Meier and Johnson, 1962). The 1982-83 surge of 
Variegated Glacier was particularly well studied (Kamb and 
others, 1985; Raymond, 1987). A surge front was observed 
on Variegated Glacier that was described as a bulge 
propagating down-glacier which was nearly coincident with 
a peak in velocity. As the front moved down-glacier, the 
peaks in height and velocity increased in amplitude and 
approached the leading edge of the wave. This resulted in a 
shock-like surge front with longitudinal compression ahead 
of the surge front and extension behind it. More recently a 
surge front was observed, via repeat image feature tracking, 
on Kunyang glacier in the Karakoram mountains of Pakistan 
(Quincey and others, 201 1 ). Over a 4 year period from 2006 
to 201 0 it was possible to track a velocity front as it grew in 
intensity and moved down-glacier, until it eventually 
diminished as the surge ended. However, the rate at which 
the surge front propagated down-glacier was not measured. 

In general, the response of Alaskan glaciers to the 
warming climate has been to retreat (Hall and others, 
2005), thin (Abalgeirsdottir and others, 1998; Arendt and 
others, 2002; Luthcke and others, 2008; Berthier and others, 
2010) and decelerate (Heid and Kaab, 2012a). Exactly how 
surging glaciers will respond to warming trends, either by 
increasing or decreasing surge frequency and magnitude, is 
unknown. A complete surge cycle, including the quiescent 
and active phases, may last from several decades to more 
than a century (Post, 1969). Therefore, it is important to 
study every surge possible, because they occur only 
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Fig. 1 . Location of Bering Glacier. Black rectangle in (a) indicates 
location of (b) within Alaska. 


intermittently, and to establish a baseline with which future 
observations may be compared to determine any possible 
changes in surge behavior. Consequently, the aim of this 
study is to observe the build-up and movement of a surge 
front on Bering Glacier as it progresses down-glacier, 
measure its speed of propagation, and compare the results 
with previous observations of the 1993-95 surge. 

METHODS 

For this study, two different feature-tracking algorithms were 
applied to sequential Landsat-7 Enhanced Thematic Mapper 
Plus (ETM+) images of the ablation zone to measure the in- 
crease in ice surface velocity through time. The first method 
used orientation images, which are produced by calculating 
the gradient in brightness values in the north-south and 
east-west directions across the satellite image, thereby 
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creating a complex image in which the east-west gradient 
is the real part and the north-south gradient is the imaginary 
part. Then the signum function (x/|x|) is applied to produce 
an image with values ranging from -1 to 1 (Fitch and others, 
2002; Haug and others, 201 0). Orientation images were used 
because they have been shown to be effective in feature 
tracking on large glaciers, specifically ice shelves in Ant- 
arctica (Haug and others, 2010), and orientation images are 
not affected by seasonal changes in solar illumination, which 
is beneficial when performing feature tracking on images 
taken during different seasons of different years, such as the 
July 2009 and September 2010 images used in this study 
(Table 1). Also, orientation images are not impacted by areas 
of uniform brightness, such as featureless snowpatches found 
in some of the winter images also used in this study (Table 1 ). 
The orientation images were inserted into COSI-Corr feature- 
tracking software (Leprince and others, 2007) to produce 
velocity maps of Bering Glacier for dates before failure of 
Landsat-7's scan-line correction mirror in May 2003. COSI- 
Corr makes use of the Fourier shift theorem which states the 
shift between two images is found in their phase difference, 
as opposed to the similarity of their intensities as with cross- 
correlation. Transformation of the images into the frequency 
domain allows isolation of the phase difference and 
measurement of the displacement between the two images 
(Shekarforoush and others, 1996). COSI-Corr was chosen 
due to its proven precision (1/10 pixel) (Scherler and others, 
2008) and because it executes faster than statistical cross- 
correlation which is performed in the spatial domain. 

The second feature-tracking program is a statistical cross- 
correlation algorithm that operates in the spatial domain 
based upon the work of Ahn and Howat (201 1). It is robust 
when given Landsat-7 ETM+ images that contain scan-line 
data voids (SLC-off images), whereas COSI-Corr is not. To 
achieve this robustness, pixels that lie within a data void are 
excluded from the cross-correlation calculations; therefore, 
only actual recorded brightness values contribute to the 
displacement measurements. Ahn and Howat (201 1 ) suggest 
the use of large reference areas with a minimum size of 
100x100 pixels to provide a sufficient number of valid 
pixels (in the presence of data voids that can be 30 pixels 
wide) to obtain accurate correlation results. This method 


Table 1 . Pairs of Landsat-7 ETM+ panchromatic imagery and feature-tracking search parameters used in this study. Search and reference 
window sizes are the same when using orientation correlation, so only one size is listed. When using statistical cross-correlation, the search 
window is larger than the reference window, so two values are given. The mean georeferencing error is given for each image pair; this value 
represents the average error in positional accuracy between the images 


Image dates 

Use in this study 

Search and reference 
window sizes 

pixels 

Mean georeferencing 
error between images 

m 

19 Apr. 2001; 22 Apr. 2002 

Orientation correlation feature tracking in lower ablation zone 

64 x 64 

4.6 

10 Sept. 2001; 29 Sept. 2002 

Orientation correlation feature tracking in upper ablation zone 

64 x 64 

2.1 

25 Apr. 2003; 1 1 Apr. 2004 

Orientation correlation feature tracking in upper and lower 
ablation zone 

64x64 

10.2 

7 Aug. 2006; 10 Aug. 2007 

Statistical feature tracking in upper and orientation correlation 
in lower ablation zone 

Upper: 128x 128, 
32 x32 

Lower: 64 x 64 

5.2 

22 Apr. 2008; 1 Apr. 2009 

Orientation correlation feature tracking in lower ablation zone 

256x256 

1.3 

30 Jul . 2009; 19 Sept. 2010 

Statistical feature tracking in upper and orientation correlation 
in lower ablation zone 

Upper: 128x 128, 
32 x32 

Lower: 256 x 256 

5.8 

19 Sept. 2010; 8 Oct. 2011 

Manual feature tracking in lower ablation zone 

Not applicable 

Not applicable 
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was used to produce velocity measurements for areas of the 
glacier affected by scan-line voids. The piedmont lobe and 
much of the ablation zone lie within an area in the satellite 
images in which the scan-line data voids disappear. So, in 
this area the combination of orientation images and COSI- 
Corr was used to produce velocity measurements, instead of 
statistical cross-correlation, after May 2003. 

Individual point displacements are assigned to the center 
of the search window from which they were derived. In the 
case of statistical cross-correlation this is reasonable 
because the displacement measurement represents how 
well pixels within the smaller reference window (centered 
within the search window) match an equally sized patch 
within the larger search window. Thus, the displacement 
represents movement of the ice encompassed within the 
reference window, taken from the first image, to a new 
location within the second image. Because the displacement 
measurement begins at the center of the windows, it is 
reasonable to assign it to this position within the first image 
of the image pair. For the case of phase correlation, such as 
COSI-Corr, the search and reference windows are the same 
size, so the dominant feature that is tracked may lie 
anywhere within the windows. It is not feasible to visually 
examine every match produced by a phase correlation 
program and determine the location of the feature matched 
within the windows; therefore, the displacement measure- 
ment is assigned to the center pixel of the search window. 

Traditionally, optical feature tracking has been performed 
on the ablation areas of alpine glaciers using summer 
images in which there are numerous ice surface features to 
track. However, initial application of the feature-tracking 
programs to summer images obtained in 2001 and 2002 
revealed the effects of shifting surface features caused when 
sediment layers emerge from the ice due to ablation (Bruhn 
and others, 2010). This happens because the exposed edge 
of the sediment layer lies within a different plane than the 
ice surface and therefore moves relative to the ice surface as 
ice melts. This phenomenon is so prevalent on the lower 
portion of the ablation zone that feature tracking using 
summer images from 2001 and 2002 yielded unusable 
results, because the ice is nearly stagnant during the 
quiescent phase and the apparent motion of the emergent 
sediment layers is greater than the down-glacier motion of 
the ice. Therefore, feature tracking was performed on 
images in which the ablation zone was snow-covered (April 
2001 and April 2002) to avoid this problem. The snow in the 
ablation zone is sufficiently deep to hide the effects of 
emergent sediment layers, while other features such as 
prominent medial moraines and crevasses have greater 
visual contrast, even when snow-covered, and are distin- 
guishable. Therefore, these moraines and crevasses can be 
reliably tracked because they dominate the local search 
window, rather than the emergent sediment layers that are 
buried by snow in winter. 

Five velocity maps were produced from the Landsat-7 
images, with time frames spanning 2001-02, 2003-04, 
2006-07, 2008-09 and 2009-10 (see Table 1 for specific 
image dates and their use). By summer 2011 the glacier 
surface was too disrupted by crevasses and contortion of 
surface features to extract accurate displacement measure- 
ments using repeat image feature tracking. Therefore, manual 
feature tracking was used to measure the displacement of 
large, obvious moraine features for 2010-11 by visually 
determining their movement from one image to the next. 
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The point displacements produced using the feature- 
tracking programs are post-processed to eliminate erroneous 
matches, remove image-to-image georeferencing error, and 
then turned into velocity rasters. Displacements on dry, 
cloud-free, snow-free land are averaged to determine the 
mean georeferencing error between sequential satellite 
images (Table 1). This error is subtracted from the on-ice 
displacements to improve their accuracy. The on-ice 
displacements are filtered using a neighborhood analysis 
routine to remove statistical outliers. Vectors that have either 
a direction or magnitude more than ±2 standard deviations 
away from their local neighborhood mean (computed from 
at least the nine nearest vectors) are deemed anomalies and 
removed. The vector field is then visually inspected and any 
remaining anomalous vectors are manually removed. The 
point data are interpolated using a linear inverse distance 
weighting scheme and then smoothed using a mean filter to 
produce a velocity raster (Fig. 2) in which the mean is 
calculated over the area contained by the search window 
used to produce the velocity field. 

The accuracy of optical feature-tracking methods depends 
upon how well the two images are co-registered to one 
another, and the precision of the algorithm used. Typical 
misalignment between two Landsat-7 ETM+ images ac- 
quired within the same World Reference System path and 
row on different days is <5 m (Lee and others, 2004). We 
find misalignments between sequential satellite images 
ranging from 1.3 to 10.2m (Table 1), with a mean (±1 std 
dev.) for all image pairs of 4.8T3.1 m, which is in 
agreement with the value give by Lee and others (2004). 
The precision of COSI-Corr is ~1 m when applied to ETM+ 
imagery, and for statistical cross-correlation it is ~9 m (Heid 
and Kaab, 2012b). Using the root-sum-of-squares method, 
the overall error in the resulting velocity maps is ±5 m when 
using COSI-Corr and ±10m when statistical cross-correl- 
ation is used. These error estimates are valid on dry land 
without deformation of surface features from one image to 
the next. Removal of the mean georeferencing error, as 
described above, will decrease these error estimates. The 
precision of the 2010-11 displacements obtained using 
manual feature tracking is estimated to be ±2 pixels (±30 m) 
due to the diffuse nature of the moraines that were matched 
and deformation of the features from one image to the next 
during the surge climax. The error in the displacement 
measurements on the glacier ice produced by the feature- 
tracking programs is difficult to quantify due to compressive 
and extensive deformation, rotation, emergent features, and 
crevassing. Statistical cross-correlation is robust against 
deformation (compression and extension) of surface features 
due to its pixel-by-pixel correlation process; however, the 
degree of its robustness has not been quantified. Phase 
correlation is more sensitive to feature deformation than 
statistical cross-correlation, but COSI-Corr incorporates a 
least-squares routine so it is also robust versus surface 
compression and extension, as well as rotation. But again, 
the degree of its robustness has not been quantified. Both 
feature-tracking algorithms are susceptible to mismatches 
from emergent features and crevassing, but these can be 
removed with the filtering routines described above. As ice 
velocity increases during the surge, it is reasonable to 
assume the degree of surface deformation will also increase. 
Therefore, we assume the amount of error in the feature- 
tracking results increases linearly with velocity, from a 
minimum of ±10m to a maximum of ±30 m, coinciding 
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Fig. 2. Velocity fields for Bering Glacier derived by repeat image feature tracking, (a) Lower ablation zone: 19 April 2001 to 22 April 2002; 
upper ablation zone: 1 0 September 2001 to 29 September 2002. (b) 25 April 2003 to 1 1 April 2004. (c) 7 August 2006 to 10 August 2007. 
(d) 22 April 2008 to 1 April 2009. (e) 30 July 2009 to 1 9 September 2010. Transect A-A' within each panel shows location of velocity profiles 
in Figure 3. Note: velocity quantization is different for each velocity field to better highlight the spatial structure within each field. 


with little deformation for slowly moving ice and with the 
error associated with visual feature tracking during the surge 
climax, respectively. Under this assumption, the error 
associated with a velocity of 1500ma _1 (roughly the peak 
velocity in 2008-09 (Fig. 3)) is ±1 7.5 m a^ 1 , or slightly more 
than ±1 pixel a -1 . At a velocity of 4390 m a -1 , the maximum 
velocity found using visual feature tracking, the error is 
±30 m a -1 , or ±2 pixels a -1 . These error estimates are ~1 %, 
or less, of their associated velocities and are smaller than the 
line thickness and size of the individual data points 
displayed in Figure 3; therefore, no error bounds are shown. 

RESULTS 

Five velocity fields were produced spanning 2001-10 
(Fig. 2). A velocity profile along transect A-A' for each of 
the five velocity fields illustrates the progression of the surge 
front down-glacier through time (Fig. 3). Beginning with the 
2001-02 profile, we interpret the small rise in velocity at 
~6 km from the confluence with Bagley Ice Valley as being 
the first observable instance of the surge front. In later years 
it is evident by looking at the peak velocities that the surge 
front steadily progresses down-glacier and the ice steadily 
accelerates year-by-year, until there is a drop in maximum 
ice velocity from ~1 .4 ±0.01 6 km a -1 in 2008-09 to 
1 .2 ± 0.01 5 km a -1 in 2009-1 0. It should be noted that the 
2008-09 velocity field was produced using winter images 
(see Table 1 ), so the velocity field does not extend up-glacier 


as far as the others due to thicker snow cover at higher 
elevations obscuring features. 

Due to large ice surface deformation and heavy 
crevassing in summer 2011, repeat image feature tracking 
was unsuccessful. However, manual tracking of medial 
moraines was still possible in the lower half of the ablation 
zone. Manual feature-tracking results, shown as individual 
points in Figure 3, reveal the piedmont lobe experienced 
nearly a threefold increase in velocity compared to the 
previous 2008-09 maximum, with velocities close to 
4.4±0.03 kma~' near the glacier terminus and decreasing 
rapidly up-glacier. 

By plotting the position of the surge front's peak velocity 
vs time it is possible to determine the rate at which the 
surge front moved down-glacier (Fig. 4). The slope of the 
line fitted through the points from September 2002 to April 
2009 gives the mean celerity of the surge front during this 
time, 4.4 ±2.0 km a -1 . This velocity is greater than the 
velocity at which the ice moves, which indicates the surge 
front moves as a kinematic wave. From April 2009 to 
September 2010 the celerity of the surge front increases to 
13.9 ±2.0 km a~\ Selecting the location of the peak of the 
surge front to track its movement down-glacier is subjective, 
with unknown errors, because the peak is not always 
obvious. Therefore, the error given for the kinematic wave 
celerity, ±2.0 km a -1 , is the standard deviation of the 
regression analysis used to determine the mean celerity 
and should be considered a minimum. 
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Fig. 3. Velocity of Bering Glacier along transect A-A' for all five velocity fields (colored lines). The surge front is seen in each profile as a step 
change in velocity, and the peak of the surge front for each profile is indicated by colored arrows. The changing location of the peak 
indicates the yearly propagation of the surge front down-glacier. The surge front increases in magnitude year-by-year from 2001 to 2009, 
until in the 2009-10 profile there is a drop in peak velocity. The surge climax is illustrated by the manually determined measurements 
(individual black squares), in which the velocity approaches 4.4 km a -1 . Error estimates for velocity are within the line thickness. Colored 
crosses on the abscissa indicate position of the glacier terminus in the second image for each image pair used to produce the associated 
velocity profile of the same color. 


DISCUSSION 
Ice dynamics of the surge 

Burgess and others (2012) found, via SAR speckle tracking 
over monthly repeat intervals, that the surge had two phases 
of acceleration surrounding a slower phase. The first 
acceleration phase, from September 2008 to February 
2009, had ice velocities up to 7md _1 (2.5 km a - '), and the 
second phase of acceleration, summer 201 1, had velocities 
up to 9md _1 (3.2 km a -1 ). These speckle-tracking results 
are of comparable magnitude to the velocities presented 
here for the April 2008— Apri I 2009 (1 .5 ± 0.03 km a' 1 ) 
and September 201 0-October 201 1 (4.4 ± 0.03 km a -1 ) time g 
frames. The slower phase described by Burgess and others, 
from January to April 2010, is within the time-span of our 
July 2009-September 201 0 velocity field and would explain ~ 
the overall decrease in ice surface velocity up-glacier of the 
surge front compared to the earlier April 2008-April 2009 
velocity field. Down-glacier of the surge front in the 2009- 
ID velocity field, the ice accelerated to ~350ma -1 , 
compared to values of 15-20ma _1 for previous years. This 
acceleration is likely due to the close proximity of the surge 
front to the glacier terminus and the associated increased 
longitudinal stress transfer down-glacier overcoming basal 
drag near the terminus. In spite of this acceleration, the 
terminus did not advance, likely due to increased calving 
into Vitus Lake. In fact, the terminus of Bering Glacier has 
retreated annually since the end of the previous surge in 
1995 (Shuchmann and others, 2010). It is only since the 
surge front reached the terminus in 201 1 that the terminus 
has advanced, roughly 2-4 km (see colored crosses in 
Fig. 4). 

The results presented here are annual measurements 
derived from satellite images acquired ~1 year apart, and 
therefore represent an average of the seasonal fluctuations in 
the surge. Conversely, Fatland and Lingle (1998) used SAR 
interferometry to measure displacements over 3 days for the 


1993-95 surge, and these more closely represent seasonal 
velocity in Bagley Ice Valley and are not directly compar- 
able to our measurements. However, visual feature tracking 



Fig. 4. Location of the peak of the surge front through time. The 
abscissa values are the date of the second image within each image 
pair used to create a velocity field; the date is given next to each 
point. The ordinate values are the distances of the maximum ice 
velocity nearest the surge front from the confluence with Bagley Ice 
Valley. The slope of the solid line fitted through the points from 2002 
to 2009 represents the average rate at which the surge front moves 
down-glacier, 4.4 km a -1 . Dotted lines represent ±2.0 km a” 1 , the 
standard deviation of the surge front speed. 
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(Roush and others, 2003) using 35 day repeat SAR images 
showed typical ice velocities of 1 0-20 m d _1 in the piedmont 
lobe between August and September 1 993. If these rates are 
extrapolated to an annual velocity (3. 7-7.3 kma -1 ), they are 
comparable in magnitude to the maximum velocity of the 
ice presented here, 4.4 ±0.03 kma~\ 

Kinematic wave 

It has been suggested that a surge front represents the 
transition of the basal hydraulic system from fast, efficient 
tunnel drainage that promotes ice movement by deforma- 
tional flow ahead of the front, to a high-pressure linked- 
cavity system behind the front that promotes flow by sliding 
(Kamb and others, 1985). Fowler (1987) describes the 
transition in terms of activation waves that move up- and 
down-glacier from a nucleation point, and he states that the 
passage of the wave indicates collapse of the tunnel drainage 
system. Fast-flowing ice behind the surge front that is 
adjacent to slow-moving ice ahead of the front creates large 
compressive stress and strain gradients across the surge front 
that induce an increase in ice thickness. Kinematic wave 
theory predicts this perturbation in ice thickness will diffuse 
over time as the wave moves down-glacier. Diffusion of the 
perturbed mass counterbalances the effects of increased ice 
velocity to stabilize the kinematic wave and produce a wave 
of constant height that moves with a constant velocity 
(Johnson, 1968). As shown above, the kinematic wave on 
Bering Glacier propagates down-glacier at an average 
velocity of 4.4 ±2.0 km a -1 between September 2002 and 
April 2009, suggesting in this instance that diffusion acts to 
stabilize the wave, as theorized by Johnson (1968). 

The kinematic wave accelerates to 1 3.9 ± 2.0 km a -1 
between April 2009 and September 201 0. This rate is derived 
from only two data points, and the straight line connecting 
them has no inherent statistical deviation, such as would be 
expected with a linear regression through a cluster of points. 
So we assign the rate the same error bounds as found for the 
wave speed from 2002 to 2009. The acceleration of the wave 
suggests a breakdown of the linear relationship between 
wave height and speed and replacement by a different 
response as it enters the piedmont lobe. Several physical 
variables change from the part of Bering Glacier contained in 
the valley to the piedmont lobe. The ice is thinnest in the 
piedmont lobe (Conway and others, 2009), the valley walls 
diminish so the ice spreads laterally to form a broad fan, the 
terminus calves into Vitus Lake, so there may be unknown 
lakewater effects, and the trend of underlying geologic 
structures changes. The orientation of geologic and topo- 
graphic structures in the mid- and upper ablation zone 
causes mountain ridges adjacent to the glacier to plunge 
beneath the ice, creating obstacles to flow. This pattern 
continues down-glacier and into the eastern half of the 
piedmont lobe. In the western half of the piedmont lobe, the 
trend in geologic structures changes to a north-south 
orientation, parallel to ice flow (Bruhn and others, 2010). 
Additionally, the existence of subglacial troughs has been 
noted in SAR imagery of the piedmont lobe (Bruhn and 
others, 2010). Exactly how all these physical attributes com- 
bine to influence kinematic wave speed is still unknown, but 
the relationship between wave height and speed has changed 
from the linear one observed up-glacier, to a new, possibly 
nonlinear relationship in the piedmont lobe. 

Between 2000 and 2003, airborne laser altimetry data 
showed there were small acceleration events in the accumu- 
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lation zone that transferred ice down-glacier to the upper 
ablation zone (Burgess and others, 2012). This formed a 
reservoir area from roughly the 30 km mark (Figs 2 and 3) 
extending up-glacier to the confluence with Bagley Ice 
Valley and eastward into Bagley Ice Valley another 5 km, 
which Burgess and others (2012) suggest acted as the trigger 
area for the first acceleration phase of the 2008-1 1 surge. In 
this discussion we refer to the location at which the surge 
began as the trigger point, or trigger area, and we refer to the 
location at which the surge front (kinematic wave) began as 
its nucleation point, or nucleation area. It is probably not a 
coincidence that our first observation of the kinematic wave 
in 2001-02 is located in the reservoir area, because Fowler 
(1987) predicts the wave will nucleate within a reservoir of 
ice that exceeds its threshold of stability. Fowler (1987) also 
predicts the surge will begin after the kinematic wave has 
propagated up- and down-glacier of the entire reservoir area 
and its hydraulic system has been activated. Our obser- 
vations seem to confirm this theory. In the current instance, 
the reservoir extends to roughly the mid-ablation zone by 
2007 (cf. fig. 5b in Burgess and others, 2012), and our 
velocity maps show the kinematic wave in this same region, 
near the 33 km mark, in 2007. The surge begins the next year 
in 2008, after the wave has passed down-glacier of the 
reservoir. Roush and others (2003) place the trigger point in 
1 993 near the 40 km mark. If it is assumed the 1 993-95 surge 
was triggered within a reservoir of ice located in the same 
region as the current surge, then the 1993-95 surge began 
down-glacier of the reservoir area as well. Roush and others 
(2003) also note the surge must have initiated up-glacier of 
the extent of a 26 March 1 993 SAR image because all the ice 
within that image was rumpled. This places the nucleation 
point for the 1 993-95 surge at least 50 km up-glacier of the 
terminus near the confluence of Bering Glacier and Bagley 
Ice Valley, and in the same general location as our first 
observation of the surge front in the 2001-02 velocity field. 

The speed of the surge front during the 1993-95 surge 
was measured using differential SAR interferometry in winter 
1992 and winter 1994 by Fatland and Lingle (1998). They 
found the surge front celerity to be up to lOOmd -1 
(36.5 kma -1 ). This rate is significantly greater than the 
estimates of 4.4 ±2.0 km a -1 and 1 3.9 ±2.0 km a -1 pre- 
sented here. Assuming the 1993-95 and 2001-10 surge 
fronts actually moved at similar rates (because of similar 
trigger and nucleation points), the difference in magnitude 
between the two surge front propagation rates may illustrate 
the difference between seasonal and annual measurements 
and could indicate a seasonal cycle to surge front propa- 
gation, faster in late winter and early spring when creep 
closure pressurizes englacial and subglacial passages, and 
slower in summer when channelization of subglacial drain- 
age pathways reduces water pressure. This type of seasonal 
pattern of surge front acceleration and deceleration was 
observed on Variegated Glacier in 1982-83 (Raymond, 
1987) in which the leading edge of the surge front moved 
down-glacier fastest during April and May, and slowest from 
July to October. Fatland and Lingle (1998) also note the 
1993-95 surge propagated up-glacier into Bagley Ice Valley 
at 200-500 md -1 (73-1 82 km a -1 ), much faster than the 
down-glacier rate. Due to limitations of optical feature 
tracking in snow-covered areas, we do not have velocity 
measurements in Bagley Ice Valley to constrain the up- 
glacier propagation of the surge and cannot make compar- 
isons with Fatland and Lingle's observations. Our obser- 
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vations show that a kinematic wave took from 2001 to 201 1 
to travel ~64 km in the ablation zone. Currently no evidence 
has been presented to show the surge front for the 1 993-95 
surge had a similar travel time, but considering the recent 
surge and the 1993-95 surge had similar trigger and 
nucleation areas, it is not unreasonable to think the two 
events had surge fronts with similar travel times. 

SUMMARY 

It has been shown that it is possible to track the build-up and 
movement of a surge front on a large temperate Alaskan 
glacier using a combination of repeat image feature-tracking 
algorithms and Landsat-7 ETM+ imagery, including images 
with scan-line voids. Analysis of the resulting velocity maps, 
spanning 2001-10, shows the surge front moved down- 
glacier in the form of a kinematic wave with an average 
velocity of 4.4 ±2.0 km a -1 between September 2002 and 
April 2009. The small variability in speed of the wave during 
this time suggests it may have been stabilized by diffusion. 
The wave then accelerated to 1 3.9 ±2.0 km a -1 between 
April 2009 and September 2010 as it entered the piedmont 
lobe. The surge appears to have climaxed in summer 201 1, 
with the ice velocity approaching 4.4 ± 0.03 km a -1 near the 
terminus. The kinematic wave is estimated to have nucleated 
near the confluence of Bering Glacier and Bagley Ice Valley 
as early as 2001. The surge began in 2008 after the 
kinematic wave moved down-glacier of an ice reservoir 
area in the mid-ablation zone, suggesting it was triggered 
there after the reservoir's basal hydraulic system was 
converted to a high-pressure linked-cavity system. 
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